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ABSTRACT 



We present results from combining a grid-based radiative transfer code with a 
Smoothed Particle Hydrodynamics code to produce a flexible system for modelling 
radiation hydrodynamics. We use a benchmark model of a circumstellar disc to de- 
termine a robust method for constructing a gridded density distribution from SPH 
particles. The benchmark disc is then used to determine the accuracy of the radiative 
transfer results. We find that the SED and the temperature distribution within the 
disc are sensitive to the representation of the disc inner edge, which depends critically 
on both the grid and SPH resolution. The code is then used to model a circumstellar 
disc around a T-Tauri star. As the disc adjusts towards equilibrium vertical motions 
in the disc are induced resulting in scale height enhancements which intercept radia- 
tion from the central star. Vertical transport of radiation enables these perturbations 
to influence the mid-plane temperature of the disc. The vertical motions decay over 
time and the disc ultimately reaches a state of simultaneous hydrostatic and radiative 
equilibrium. 
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1 INTRODUCTION 

Hydrodynamics calculations are used in investigations of a 
variety of problems associated with star formation, such as 
the collapse of protostellar cores (e.g. Goodwin et al.|2004| ), 
the fragmentation of protostellar discs (e.g. |Boss|1997| ), and 
on the largest scale the gravo-turbulent collapse of a molecu- 
lar cloud to form stars (e.g. Bate et al. 2003 Bate &; Bonnell 
[2005 ). 

Sophisticated as such calculations appear, the under- 
lying physics regarding the thermal properties of the gas 
is much simplified. One of the usual approximations is the 
use of an equation of state (eos) in determining the pres- 
sure of the gas. An isothermal eos is often a good approx- 
imation in the most rarified regions, where the material is 
able to radiate, and thus cool, freely. The gas may be com- 
pressed adiabatically in the highest density regions, and a 
barotropic eos is sometimes adopted, with the poly tropic 
exponent switching at a critical density in an attempt to re- 
produce the results of one-dimensional collapse calculations 
that incorporate radiative-transfer |Masunag a et al.p 998). 
The crucial point here is that the energy dumped into the 
cluster from accretion onto, and gravitation contraction of, 
the newly formed stars is neglected, and the radiation emit- 
ted by compressed gas (the L PdV term) is assumed to es- 
cape without further interaction. 



The solution is to include radiation-transfer self- 
consistently and perform a radiation- hydrodynamics (RHD) 
calculation. The simplest way to do this is to assume that the 
energy is transported by radiative diffusion of photons. The 
diffusion approximation is a good one at high optical depth, 
but breaks down in optically thin regions since the diffu- 
sion speed may exceed the speed of light as the mean-free- 
path tends to infinity; an essentially arbitrary flux limiter is 
adopted in this regime. Furthermore the approximation is 
often implemented as a grey method, typically using Rosse- 
land mean opacities tabulated as a function of temperature. 
Finally the diffusion approximation cannot treat the effects 
of shadowing, as the radiation field may diffuse around op- 
tically thick obstacles. The crucial advantage of the diffu- 
sion approximation though is of course speed: the radiation 
transport step is straightforward to implement and is not 
dominant in the RHD calculation. The first cluster collapse 
and disc fragmentation calculations using flux-limited diffu- 
sion have just started appearing in the literature flBoss|2008| 
[Bate]|2009bl |Offner et aTl|2009| |Stamatellos fe Whitworth| 
2009). The fundamental result of including a more complete 
thermal treatment is broadly in line with prior expectations 
i.e. the heating of the gas inhibits the fragmentation process. 

Meanwhile dedicated radiation transfer (RT) codes have 
progressed far beyond the diffusion approximation and now 
take into account a wide array of processes such as poly- 
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chromatic radiation transfer incorporating the radiation's 
polarization state, multiple scattering, Mie phase matrices 
etc. These sophisticated RT codes have been employed to 
model co res (e.g. |Stamatellos et al. 2005), discs (reviewed 
by Dullemond et al. 2007) and even entire clusters ( |Kuro-| 
sawa et al.|2004 ). The principle barrier to using these codes 
in RHD is that of of computation time, but the increase in 
CPU power and the adoption of parallelization mean that it 
is now feasible to conduct a full treatment of RT along with 
the hydrodynamics, thereby dropping the diffusion approx- 
imation. 

This paper presents a method of coupling the TORUS 
radiative transfer code (Harries 2000|) with the smoothed 



particle hydrodynamics (SPH) code of Bate (2009a). TORUS 



performs radiative transfer calculations using the Monte- 
Carlo method of |Lucy| ( |1999| ) with a diffusion approximation 
in regions of high optical depth. The radiative transfer cal- 
culations are performed on a grid which uses adaptive mesh 
refinement (AMR) to allow variable resolution. The SPH 
method is well suited to dealing with the large range in size 
scales encountered during star formation processes, hence 
an SPH code coupled to a Monte-Carlo radiative transfer 
code is a very flexible method for performing star formation 
calculations incorporating radiative feedback. 

Accurately coupling the radiative transfer component 
to the SPH component requires an effective transformation 
from the Lagrangian particle description of SPH to the Eu- 
lerian grid description of an AMR grid. The procedure used 
for this transformation is tested using a circumstellar disc 
benchmark (see Sections [2] and [3| . The temperature distri- 
bution and spectral energy distribution of the benchmark 
disc are examined in Section [4] in order to investigate how 
well the radiative transfer calculation works with the grid- 
ded density field and to understand the effects of grid and 
SPH spatial resolution. Finally the hybrid code is tested un- 
der realistic conditions via a simulation of a circumstellar 
disc around a T Tauri star (Section [5]). 



2 CIRCUMSTELLAR DISC BENCHMARK 

The accuracy of the SPH method has previously been ex- 
amined in detail (e.g. [Price (2008) and references therein) so 
we do not attempt to benchmark the SPH component of the 
code here. Likewise the TORUS radiative transfer code has 
also been benchmarked elsewhere ( |Pinte et al.||2009| ). Con- 
sequently we can be confident in the accuracy of the SPH 
and radiative transfer codes individually and we instead fo- 
cus on the accuracy of the mechanism by which they are 
coupled. The benchmark case presented here is therefore a 
'static' benchmark in which we represent an analytical den- 
sity distribution with SPH particles and test whether the 
radiative transfer code can calculate accurate results using 
this SPH density representation. 

In the absence of an analytical test case, benchmark- 
ing of codes used for calculating radiative transfer in cir- 
cumstellar discs is carried out by code comparison studies 
( |Pascucci et~aLl|2004| |Pinte et al.1|2QQ9| ). We base our SPH 
disc benchmark on the axisymmetric circumstellar disc of 
( 2004] ). 



Pascucci et al. (20041). Representing a structure with rota- 
tional symmetry using a 3D geometry allows us to validate 
the operation of our method against well tested results from 



other codes while fully exercising the 3D capabilities of our 
own system. The geometry is a stringent test of the grid- 
ding method, since it contains a large range in linear scales 
that need to be resolved, and sharp density gradients (in 
particular the disc inner edge). 

The density of the disc benchmark described by |Pas-| 
cucci et al. ( 2004 ) is given by 



p (r, z) p 



exp 




(i) 



where r is the distance from the central star in the mid-plane 
and z is the height from the mid-plane. The expression h (r) 
is given by 



h(r) 



Zd 



(2) 



and po, Zd and Td are constant parameters. The definition of 

J200l 



h (r) is the definition used by |Pascucci et al.| ( |2004| which 
differs from the standard scale height definition (used later 
in this paper) by a factor of f. We base our benchmark 
model on the most optically thick case presented by Pascucci 
et al. (| 2004| which has a mid-plane optical depth of r — 100 
550 nm and a mass of 0.011 Mq. This requires rd — 
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500 au, z d = 0.25r d and p = 8.1614 x 10 _1 *gcm _;3 , with 
the disc truncated at an outer radius of 1000 au and an inner 
radius of 1 au. The disc is irradiated by a solar-type star with 
a black body spectrum. 

The benchmark disc is implemented using an ensemble 
of equal mass SPH particles. If all the particles have the 
same mass then the number density is proportional to the 
mass density, hence the SPH particles sample a probabil- 
ity density distribution such that the probability of finding 
a particle in a given volume is proportional to the mass in 



that volume as a fraction of the total mass ( Gingold & Mon- 
aghan 1977). Hence the analytical density function can be 



converted into a probability function which describes the 
probability of finding a particle in a given volume. An en- 
semble of particle positions can then be calculated by ran- 
domly sampling the probability density functions for r and 
z and assigning a uniformly random azimuthal angle 0. The 
particle's density is set to the value given by Eqn[l] based 
on the particle's position. 

Each particle is assigned a smoothing length h smoo th 
given by 



h s 



1.2 



\ Ppart J 



1/3 



(3) 



where m pa rt is the particle mass and p par t is the particle 
density, according to the method of Price &; Bate (2007). 
The smoothing length is used when reconstructing proper- 
ties represented by the ensemble of particles, for example 
the density distribution on the radiative transfer grid (see 
Section [3| . 



3 GENERATION OF AN ADAPTIVE GRID 
FROM SPH PARTICLES 

An important step in setting up grid-based radiative transfer 
using a density distribution derived from SPH particles is de- 
termining an effective transformation from the Lagrangian 
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particle description to the Eulerian grid description. The 
representation on grid should preserve important properties 
of the SPH model (e.g. total mass) but also needs to take 
account of the possibility that a radiative transfer calcula- 
tion requires the highest resolution in different regions to 
a hydrodynamics calculation (e.g. to resolve opacity gradi- 
ents). This section determines a grid construction method 
which allows a good degree of control over the number of 
cells in the grid, and the location of the highest resolution 
regions, while maintaining an accurate representation of the 
total disc mass. 



3.1 Grid generation method 

TORUS uses an adaptive grid based on the octree method. In 
a 3D geometry the initial grid comprises 8 cells (one octal) 
with 2 cells in each dimension. An algorithm is applied to 
determine whether to split each cell into a further 8 cells, 
similar to the method of Kurosawa & Hillier (2001). The 



properties of the adaptive grid are determined by a combi- 
nation of the algorithm used to decide when to split a cell 
and the method used to assign density values to cells. 

The density in a given grid cell is calculated using an 
exponential kernel smooth. If the sum of the kernel weights 
is greater than 0.3 the density is normalised by the sum of 
the kernel weights, to ensure a smooth density distribution 
within the interior of the disc. If the sum of the weights is 
less than 0.3 then the density is not normalised by the sum 
of the weights, in order to avoid numerical effects at the 
free surface, as described by Price (2007). For more details 
of the algorithm used to generate the AMR grid from SPH 



particles see Rundle & Harries ( 2009 ) 



3.2 Cartesian grid generation tests 

Two grid splitting conditions were tested. The first method 
splits a grid cell if the mass within the cell exceeds a given 
limit. For this condition the resolution is highest where the 
density is highest and the grid cell size is analogous to the 
SPH particle smoothing length. The second condition de- 
cides whether to split the cell based on the fractional den- 
sity difference between the most dense and least dense SPH 
particles contained within the cell. The quantity 



/split - — — (4J 

Pmax ~~r Pmin 

is calculated, where p max is the highest SPH particle den- 
sity and p m in is the lowest SPH particle density, and the 
cell is split if a specified value of /split is exceeded. The first 
condition gives a resolution like that of the original SPH 
representation whereas the second condition should better 
represent gradients in density which are important for ra- 
diative transfer calculations. 

To illustrate the effects of the different splitting meth- 
ods two example grids are shown in Fig. ^ A grid con- 
structed from 10 7 particle s, usin g a mass per cell limit of 



5 x 10 g, is plotted in Fig. 1(a) This plot is a slice perpen- 



dicular to the disc mid-plane and shows increased resolution 
towards the disc centre and mid-plane where the density is 
higher. A grid with / sp Ht =0.1, also using 10 7 particles, is 



also increases towards the centre of the disc. When the den- 
sity split condition is used there is a noticeable asymmetry 
to the grid resolution, due to statistical fluctuations in the 
SPH particle distribution. The effect is seen more readily 
with the density split condition because the grid resolution 
is being enhanced in lower density regions, where SPH par- 
ticles sample the density distribution less well than in high 
density regions. 

For each splitting condition a number of AMR grids 
were generated using different values of the mass per cell 
limit or / sp iit- This was repeated for discs represented by 
10 5 , 10 6 and 10 7 SPH particles. Figure ||] plots the num- 
ber of octals generated as a function of mass per cell or 
/split (solid line), and the percentage error in the disc mass 
(dashed line). The percentage error in disc mass is calcu- 
lated by comparing the total mass on the AMR grid with 
the known mass of 0.011 Mq based on the analytical form of 
the density distribution, and a positive value indicates that 
mass on the grid is too large. 

The mass per cell condition allows a wide range in the 
total number of octals while maintaining a total mass which 
is correct to within a few percent. The number of octals as 
a function of mass per cell limit is consistent between runs 
with different numbers of SPH particles, apart from the case 
with 10 5 particles and a mass limit of 10 26 g. In this case the 
mass per cell limit is smaller than the SPH particle mass and 
the grid is over sampling the SPH resolution. In the other 
cases the mass per cell limit is greater than the SPH particle 
mass. For mass per cell limits of 10 28 g and less the total 
mass error varies only slowly as a function of mass per cell 
and is accurate to within a few percent. There is a positive 
bias in the total mass but using more SPH particles results 
in a more accurate total mass. 

If a density contrast condition is used then the total 
number of octals depends on the number of particles, unlike 
in the mass per cell case. The general form of the number 
of octals as a function of density contrast is the same in 
all three cases but the normalisation varies substantially. In 
order to ensure a total mass accurate to with 1 per cent 
the density contrast condition would require at least 10 7 
particles and a density contrast threshold of no more than 
0.3. 

The density contrast condition is effective in adding 
extra resolution in regions of high density contrast which 
are likely to be important in radiative transfer calculations. 
However this condition has a less robust total mass represen- 
tation than the mass per cell limit. By using a combination 
of these two conditions it is possible to add extra resolution 
in high contrast regions and still maintain an accurate total 
mass. 



3.3 Cylindrical polar grid generation 

For some geometries it may be appropriate to use a cylin- 
drical polar co-ordinate system. A mass per cell splitting 
condition can still be used but needs to be modified to take 
into account that for fixed values of dr and dz the cell vol- 
ume increases as cylindrical polar radius increases. Consider 
a Cartesian grid with cubic volume elements and let 



plotted in Fig. 1(b) The resolution increases towards the 
edge of the disc, where there are large density gradients, and 



dlc&rt — dx — dy — dz 
so that 



(5) 
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(a) AMR grid with a mass per cell condition of 5 x 10 26 g. (b) AMR grid with a density condition of /split =0.1. 

Figure 1. Example AMR grids for a mass per cell splitting condition of 5 x 10 26 g (left) and a density contrast splitting condition of 
/split =0.1 (right) constructed from 10 7 SPH particles. The grids are slices perpendicular to the disc mid-plane. 




(a) Mass per cell condition (b) Density contrast condition 

Figure 2. Number of octals (solid line) and percentage error in total disc mass (dashed line) for the benchmark disc represented on 
Cartesian AMR grids. The grids were generated using either a mass per cell limit (left) or density contrast limit (right) and different 
numbers of particles (10 5 , 10 6 or 10 7 ). Results using 10 5 particles are plotted with squares, 10 6 particles with circles and 10 7 particles 
with triangles. The number of octals for a mass per cell condition with 10 7 particles is not plotted as it is indistinguishable from the case 
with 10 6 particles. 



dl c 



dM c; 



1/3 



(6) 



where dM cart is the mass of a cell in the Cartesian grid and 
Pceii is the density of the cell. In a cylindrical polar grid the 
mass of a cell is given by 



dM cy \ — yOcell 



27irdrdz 



(7) 



if there are n az uniform azimuthal cells. If the resolution in 
r and z is the same then we can define 

dl cy i — dr — dz (8) 

which is related to the cell mass and density by 

1/2 



\ 27rrp ce ii J 



(9) 
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If we require the resolution in cylindrical polar r and z to 
match the resolution in Cartesian x, y and z then dl cy \ = 
dlcart- From eqn[6]and[9] 

dMcyl= ^ pl /3 dM 2/3 (10) 

which describes the required conversion from a mass per 
cell limit suitable for a Cartesian grid (dM cart ) to an equiv- 
alent mass per cell limit suitable for a cylindrical polar grid 
(dMcyi). 

The grid generation tests were run using a cylindrical 
polar geometry with the modified mass per cell limit de- 
scribed above and an unmodified density contrast condition. 
The results are plotted in Fig. [3] Using a cylindrical polar 
geometry results in a significantly reduced number of oc- 
tals in the AMR grid but with larger mass errors. Provided 
the mass per cell limit is no more than 10 27 g then the total 
mass is well represented with significantly fewer cells than in 
a Cartesian grid. The density contrast condition alone does 
not represent the total mass to within 3% accuracy even 
with / S pnt =0.1. Although this condition successfully adds 
enhanced resolution in regions of higher density contrast it 
needs to be combined with a mass per cell limit to ensure a 
robust total mass representation. 



4 COMPARISON OF TEMPERATURE 
DISTRIBUTION AND SEDS WITH 
BENCHMARK RESULTS 

Having investigated the procedure for setting up the radia- 
tive transfer grid we now proceed to studying the accuracy 
of the results from the radiative transfer calculation. This 
was carried out by comparing temperature distributions and 
spectral energy distributions (SEDs) from the benchmark 
disc. Temperature distribution errors provide a clear pic- 
ture of where errors in the radiative transfer calculation are 
located whereas SEDs provide a direct link to an observable 
quantity and are a stringent test of the accuracy of the ra- 
diative transfer calculation. The key problem here is that 
the resolution necessary for adequately describing the fluid 
flow is not necessarily the same as required for accurately 
modelling the radiation field. 

The disc was initially modelled using 10 5 SPH particles 
with an AMR grid constructed from a mass per cell limit 
of 5 x 10 26 g. This number of SPH particles would normally 
be considered as sufficient to adequately represent a disc 
of this mass in a purely hydrodynamical calculation. The 
mass per cell limit was chosen in preference to the density 
contrast limit to give the highest resolution near the central 
region of the disc, which we expect to be most important for 
determining the temperature of the disc at the inner edge 
and the SED. 

In Fig. [4] we plot the fractional error in the temperature 
distribution at radiative equilibrium for the benchmark disc. 
TORUS calculates the temperature distribution using a high 
resolution 2D geometry, in order to determine the correct 
temperature as a function of cylindrical polar r and z co- 
ordinates. This is then compared to the temperature from 
the SPH disc for each SPH particle. Positive values of the 
fractional error indicate that the SPH disc is hotter than the 
expected value. The temperatures show deviations of —30% 
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Figure 4. Fractional errors in the temperature distribution of a 
10 5 particle disc. The radiative transfer grid is a Cartesian grid 
generated using a mass per cell limit of 5 x 10 26 g. Axis units are 
au. 



to +50% from the benchmark values, a deviation that is 
purely driven by the resolution of the density field derived 
from the SPH particles, since a 'TORUS only' calculation, 
where the AMR grid is refined on optical depth criteria and 
the cell densities are found from the analytical description 
of the density structure (Eqn. [T]), matches the benchmark 
temperatures to within 2%. The systematic differences in 
the temperature distribution are due to an increased verti- 
cal transport of radiation from the surface layers of the disc 
to the mid-plane due to a lack of spatial resolution, leading 
to errors in the calculation of the specific intensities (partic- 
ularly in the inner disc region) . These systematic differences 
are more readily apparent in the emergent SEDs, which are 
a more sensitive probe of the accuracy of the RT calculation 
than the temperature distribution. 

SEDs were calculated using the TORUS code for viewing 
angles of 12.5 and 77.5 degrees (where an inclination angle 
of zero degrees indicates that the disc is viewed face-on) and 
the results are plotted in Fig. [5] (dashed line). Also plotted 
in Fig. [5] are the benchmark results of |Pascucci et al.| (|2004) 
(solid line). Both SEDs show significant departures from the 
benchmark results; there is excess emission between 10-100 
/im and at a viewing angle of 77.5 degrees there is also a 
significant reduction in flux below 1 /im. Both these discrep- 
ancies are attributable to a lack of resolution within a few au 
of the central object. The density distribution in the central 
region of this disc is plotted in Fig. [6] which shows that the 
central 1 au gap is not represented. There are two effects in 
operation; firstly the grid resolution in this region is larger 
than 1 au and secondly the smoothing lengths of the parti- 
cles are larger than 1 au. In order to correctly represent the 
central gap there needs to be higher grid resolution around 
the central source and more SPH particles are required so 
that the smoothing lengths are smaller. 
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(a) Mass per cell condition 



(b) Density contrast condition 



Figure 3. Number of octals (solid line) and percentage error in total disc mass (dashed line) for the benchmark disc represented on 
cylindrical polar AMR grids. The grids were generated using either a mass per cell limit (left) or density contrast limit (right) and 
different numbers of particles (10 5 , 10 6 or 10 7 ). Results using 10 5 particles are plotted with squares, 10 6 particles with circles and 10 7 
particles with triangles. The number of octals for a mass per cell condition with 10 7 particles is not plotted as it is indistinguishable 
from the case with 10 6 particles. 
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Figure 5. SEDs from the benchmark disc represented by 10 5 particles and a Cartesian radiative transfer grid generated using a mass 
per cell limit of 5 x 10 26 g. SEDs are shown for inclination angles of 12.5 degrees (left) and 77.5 degrees (right). The benchmark SED 
is plotted as a solid line and the TORUS SED is plotted as a dashed line. The corresponding fractional temperature errors are plotted in 
Figg] 



4.1 Grid resolution 

Decreasing the mass per cell limit to achieve the required res- 
olution around the central source would result in too many 
cells being generated in other regions of the disc. In order to 
achieve a grid resolution of less than 1 au around the source 
it was necessary to specify additional criteria for refining 
the grid. Three cubes of increased resolution were added 
to the AMR grid with sizes 1.5 x 10 15 cm, 3 x 10 14 cm and 
6 x 10 13 cm. Within each cube it was required that the cell 
size be no larger than 0.1 of the cube dimension. In order to 
confirm that the grid modifications provide sufficient reso- 
lution to generate an accurate SED three calculations were 
performed with the central part of the disc forced to the cor- 



rect density values from the analytical density description. 
A spherical region with radius 1.0 x 10 15 cm, 5.0 x 10 14 cm 
or 2.5 x 10 14 cm was over written with the benchmark disc 
density after the initial grid generation was performed. The 
resulting SEDs are plotted in Fig. [7] The solid line shows 
the benchmark result and the other lines show SEDs for 
forcing regions of radius 1.0 x 10 15 cm (long dashed line), 
5.0 x 10 14 cm (short dashed line) and 2.5 x 10 14 cm (dot- 
ted line). With the correct density forced out to a radius 
of 1.0 x 10 15 cm the SED is close to the benchmark result 
for both viewing angles, however if the forced region is re- 
duced in size there is a significant discrepancy relative to 
the benchmark SED. This confirms that the grid now has 
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Figure 7. SEDs from the benchmark disc represented by 10 5 particles and a Cartesian radiative transfer grid generated using a mass 
per cell limit of 5 x 10 26 g. Extra resolution has been added around the central source and the TORUS density distribution is forced to 
the correct value within a radius of 1.0 x 10 15 cm (long dashed line), 5.0 x 10 14 cm (short dashed line) and 2.5 x 10 14 cm (dotted line). 
SEDs are shown for inclination angles of 12.5 degrees (left) and 77.5 degrees (right) with the benchmark result plotted as a solid line. 
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Figure 6. Density (gem -3 ) in the central region of the bench- 
mark disc represented by 10 5 particles and a Cartesian radiative 
transfer grid generated using a mass per cell limit of 5 x 10 26 g. 
The central 1 au gap is not represented and the effects on the 
temperature distribution and emergent SEDs are shown in Fig. ^] 
and Fig. [^respectively. 

enough resolution to generate an accurate SED but also in- 
dicates that using 10 5 particles does not give a sufficiently 
accurate density distribution within the central 1.0 x 10 15 cm 
of the disc. 



4.2 Number of particles 

Further calculations were performed using different numbers 
of SPH particles and the enhanced grid resolution described 
above in order to determine the impact of the number of 
SPH particles on the accuracy of the radiative transfer cal- 
culation. The density distribution in the central region of 



the disc is plotted in Fig. [8] in order to show how the repre- 
sentation improves as the number of particles is increased. 
The effective resolution of the AMR mesh is dependent on 
two factors; the probability of having sufficient particles to 
adequately sample the small volume of the disc inner-edge, 
and the particle smoothing length (which determines the ef- 
fective spatial resolution of the SPH density representation). 
With 10 5 and 10 6 particles a central gap is present but is too 
large. With 10 7 particles there is more material close in to 
the source but the gap itself is not represented (because the 
smoothing length of the particles near the disc inner-edge is 
too large) and the source is obscured. Only with 10 8 parti- 
cles is a gap of approximately the correct size present. The 
representation of the central 1 au gap is significantly affected 
by the number of particles used to represent the disc, with 
a large number of particles required to achieve an accurate 
density distribution in this spatially small region. 

Figure [9] shows fractional errors in the temperature dis- 
tribution for different numbers of particles. When using 10 5 
particles to represent the disc there is a too much heating in 
the disc mid plane and insufficient heating further out of the 
mid plane. The effect is reduced as the number of particles 
is increased and reaches a very low level with 10 8 particles, 
at which stage the temperature distribution agrees with the 
benchmark to within ^10%. 

To allow a more quantitative comparison the mid-plane 
temperature and the benchmark mid-plane temperature are 
plotted in Fig. ^] for all particles within 0.01 scale heights 
of the mid-plane. The fractional temperature error for each 
particle is also plotted where a positive error indicates a 
temperature hotter than the benchmark value. 

SEDs are plotted in Fig. for a viewing angle of 12.5 
degrees and Fig. [12] for a viewing angle of 77.5 degrees. At 
the shortest wavelengths the 12.5 degrees inclination angle 
SED is dominated by the photosphere of the central star. 
This region of the SED is well represented with 10 5 and 
10 8 particles but is of course less well represented with 10 6 
or 10 7 particles where the central star is obscured. At the 
longest wavelengths the disc is optically thin and the SEDs 
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Figure 8. Density (g cm 3 ) in the central region of the disc as represented on AMR grids derived from different numbers of SPH part icles. 
The AMR grid is generated using a mass per cell limit of 5 x 10 26 g with extra resolution at the origin as described in Section |4~l| 



will be well modelled provided the total disc mass is ac- 
curately represented. In all four cases the total disc mass 
is well represented (see Section [3| and the long wavelength 
SED matches the benchmark results. At intermediate wave- 
lengths the SED is sensitive to the details of the disc repre- 
sentation and it is clear that accurately modelling the SED 
from this disc is challenging, even with 10 8 particles. At an 
inclination angle of 77.5 degrees the central star is overly ob- 
scured by disc material. If the spatial resolution at the inner 
edge is low then the inner edge will be too thick (in the z 
direction) and too close to the central star. As a result there 
will be overly dense material in the observer's line of sight. 



An SED with insufficient emission from the central star at 
this inclination angle is indicative that the inner edge is not 
well resolved. 

The results of our benchmark tests have shown that 
the temperature distribution within the disc can be calcu- 
lated with a typical accuracy of better than 20 per cent 
with only 10 6 particles (see Fig. [9] and Fig. 10), which offers 
a significant improvement over the use temperatures derived 
from an eos. Although relatively few particles are required 
to give a temperature distribution sufficiently accurate for 
hydrodynamic calculations, obtaining an accurate SED from 
the model disc would require many more particles, as there 
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Figure 9. Fractional errors in the temperature distribution of discs represented by 10 5 , 10 6 , 10 7 and 10 8 parti cles. The AMR grid is 
generated using a mass per cell limit of 5 x 10 26 g with extra resolution at the origin as described in Section [XT] The axis units are au. 
The corresponding density distributions in the central part of the disc are plotted in Fig. [8] 



are significant discrepancie s rel ative to the benchmark even 
with 10 8 particles (see Fig. 11 and Fig. 12). 



are higher than that of the Pascucci benchmark, providing 
a more demanding challenge both in terms of computation 
time and required numerical accuracy. 



5 COUPLED SPH AND MONTE-CARLO 

RADIATIVE TRANSFER SIMULATION OF 
A CIRCUMSTELLAR DISC 

The above analysis indicates that the hybrid code is capa- 
ble of successfully passing the density distribution from the 
hydrodynamics step to the radiative-transfer code, and that 
temperatures sufficiently accurate for hydrodynamic calcu- 
lations (under the assumption of radiative equilibrium) can 
be passed back to the SPH code in order to update the 
gas pressure and perform the next hydrodynamical step. 
In order to test the feasibility of performing a radiation- 
hydrodynamics calculation using the hybrid code we adopt 
a 'realistic' test problem, that of a circumstellar disc around 
a classical T Tauri star. The mid-plane optical depths here 



5.1 Description of model 

A model system was constructed similar to the classical T 



Tauri star AA Tau. O'Sullivan et al. (2005) determine the 



properties of A A Tau using Monte- Carlo modelling of the 
observed SED, and we use these derived properties as a basis 
for constructing our model system. The disc of 0.02 M© is 
initialised with an outer radius of 150 au, and an inner radius 
of 1.0 au. The inner radius of our model disc is enlarged, 
compared to the inner radius determined by O'Sul livan et al.| 
(12005} , as resolving a smaller inner edge in a disc of this size 
would require many more SPH particles. The inner edge of 
our model disc will be cooler than the inner edge of a disc 
with a smaller radius but the important physical processes 
are captured (i.e. there is an optically thick inner edge which 
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can puff up and shadow the disc causing a hydrodynamic 
response.) The central star has a mass of mass 1 Mq and is 
represented by a sink particle which accretes gas particles 
within an accretion radius of 1.0 au. The disc itself is initially 
represented using 10 6 particles of equal mass, although some 
of these particles are accreted by the central sink particle as 
the simulation progresses. The central star has a radius of 
2R , a temperature of 4000 K and a luminosity of 0.9 L . 
The source spectrum taken to be a black body. 

The disc is initially evolved using only the SPH part of 
the code in order to bring the disc into dynamical equilib- 
rium. The SPH implementation is that described by |Price~fe| 
Bate ( 2007 ) which uses the variable smoothing length formu- 



lation given in eqn[3| Density and smoothing length are cal- 
culated iteratively according to the method of Pric e &; Mon-| 
aghan (20071) with approximately 60 neighbours on average. 



An artificial viscosity term is calculated using the scheme 
described by Monaghan (1992) with parameters a = 1 and 
[3 — 2. During this initial phase the internal energies of 
the SPH particles are determined using an equation of state described by Nayakshin et al. 2009). Calculating temper- 



which is isothermal perpendicular to the plane of the disc 
and has a T oc 1/r radial dependence. The disc is evolved 
for approximately 4000 years of simulated time in order to 
allow transient effects from the initialisation to dissipate. As 
the disc settles towards an equilibrium state there are verti- 
cal motions which cause fluctuations in the envelope of the 
disc, with the inner regions of the disc reaching equilibrium 
more rapidly than the outer regions. 

The disc is then evolved using the combined SPH and 
radiative transfer code. In this configuration the internal 
energies of the SPH particles are set using the tempera- 
tures calculated by the radiative transfer code. The SPH 
part of the code does not change the particle temperatures 
and they retain the temperature assigned by the radiative 
transfer code until the next radiation time step. The tem- 
peratures calculated by the radiative transfer code are equi- 
librium temperatures, as it is assumed that the time scale 
for radiative equilibrium is short compared to the dynami- 
cal time scale of the disc (i.e. the "prompt escape" regime 



Modelling circumstellar discs with 3D radiation hydrodynamics 11 




10 100 
Wavelength (jim) 




10 100 
Wavelength foxm) 



(a) 10 5 particles 



(b) 10 6 particles 




10 100 1000 

Wavelength ((im) 




10 100 1000 

Wavelength (|im) 



(c) 10 7 particles 



(d) 10 8 particles 



Figure 11. SEDs at a 12.5 degree inclination angle for discs represented by 10 5 , 10 6 , 10 7 and 10 8 particles. The benchmark result is 
plotted as a solid line and the TORUS SEDs are plotted as dashed lines. The corresponding density distributions in the central part of 
the disc are plotted in Fig. [8] and fractional temperature errors are plotted in Fig. [9] 



atures from the radiative transfer scheme alone is only a 
valid approximation provided heating from stellar radiation 
dominates other sources of heating (e.g. viscous or shock 
heating) . The addition of a time dependent radiative trans- 
fer scheme would allow for the inclusion of additional (non- 
stellar) sources of heating, such as adiabatic compression, 
but is beyond the scope of this work. 

The frequency with which the radiative transfer is cal- 
culated is determined by the sound crossing time in the inner 
regions of the disc as this is the time scale over which the 
scale height in this region will vary and effect the transfer of 
radiation into the outer parts of the disc. The vertical sound 
crossing time at a radius of 10 au is typically approximately 
2 years, so a radiative transfer calculation is performed after 
every four hydrodynamics time steps (equivalent to every 2.4 
years) . 

A cylindrical polar grid, with adaptive r and z cells was 
used for the radiative transfer step. The cells were split in 
the azimuthal direction to give a cell spacing of 22.5 de- 



grees. This azimuthal splitting allows for the existence of 
non-axisymmetric instabilities (such as radiatively driven 
warps) which may occur due to either physical or numerical 
effects, while ensuring the calculation remains tractable. The 
AMR grid was constructed using a maximum mass per cell of 
3.0 x 10 25 g and a density splitting condition of / sp iit = 0.1. 
The density in a cylindrical region around the source was 
reduced to a very low value to ensure that the disc did not 
completely obscure the source. The radius of this region is 
^Rgap = 1.0 au which corresponds to the accretion radius of 
the SPH sink particle. Two extra levels of grid refinement 
were added in a box around the central source to ensure that 
the grid had enough resolution to allow the gap around the 
source to be resolved. The first box of extra refinement has 
size 1 x 10 14 cm with cells no larger than 1 x 10 13 cm and 
the second box has size 4 x 10 13 cm with cells no larger than 
4 x 10 12 cm. 

Monte- Carlo radiative transfer calculations can be com- 
putationally demanding so the calculation was parallelised 
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Figure 12. SEDs at a 77.5 degree inclination angle for discs represented by 10 5 , 10 6 , 10 7 and 10 8 particles. The benchmark result is 
plotted as a solid line and the TORUS SEDs are plotted as dashed lines. The corresponding density distributions in the central part of 
the disc are plotted in Fig. [8] and fractional temperature errors are plotted in Fig. [9] 



using MPI (Message Passing Interface) libraries. Parallelis- 
ing code in this way allows for parallel computation in both 
the radiation and hydrodynamics steps and can be applied 
to a shared memory or a distributed memory computer ar- 
chitecture. The Monte-Carlo calculation dominates the run 
time of the model but fortunately is very amenable to par- 
allelisation, due to the low communication overhead. The 
calculations were run on an SGI Altix ICE system using 
16 compute nodes, each node comprising two quad core 
2.83 GHz Intel Xeon processors. A total of 128 MPI pro- 
cesses were used (one per core) . Evolving the disc using both 
the radiation and hydrodynamics codes took a total wall 
time of 560 hours of which 440 hours was spent in radiation 
calculations. The radiation calculations take approximately 
four times as long as the hydrodynamics calculations so the 
combined code is significantly slower than the SPH code 
alone. However the calculation is still feasible and we expect 
the combined code only to be used in applications where 
the increased accuracy of the Monte-Carlo radiation calcu- 



lation, over faster alternatives such as flux-limited diffusion, 
is of importance. 



5.2 Results 

Plots of the disc internal energy are shown in Fig. ^3] at six 
different times during the simulation. As these plots are in 
cylindrical polar co-ordinates, with different scales on the r 
and z axes, they have been shown as particle plots, rather 
than kernel integrated plots. 



Figure 13(a) shows the ini- 
tial state before the radiative transfer is switched on when 
there is a vertically isothermal equation of state. The second 



frame (Fig. 13(b) is shortly after the radiative transfer has 



been switched on. The inner disc has collapsed and there is 
heating seen in parts of the disc exposed to radiation from 
the central star. In Fig. 13(c) a scale height enhancement is 



present in the disc at a radius of approximately 20 au which 
shadows the region immediately behind causing it to cool. 



Figure 13(d) shows several such scale height enhancements 
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Figure 13. Internal energy distribution in the model disc at six different times as the disc adjusts towards equilibrium. Figure [13 (a) | 
shows the initial vertically isothermal state and Fig. |13(f)1 shows the steady state in hydrostatic and radiative equilibrium. Intermediate 
frames show the effect of scale height fluctuations on the internal energy within the disc. Internal energy is plotted for each particle in a 
cylindrical polar co-ordinate system and the axis units are au. 



and in each case there is cooling behind the enhancement 
as the scale height fluctuations are able to affect the mid- 
plane temperature due to vertical transport of radiation. In 



Fig. 13(e) the inner region of the disc, within approximately 



only in the outer disc, implying an equilibrium time of or- 
der 1600 years at 75 au. As the equilibrium time scales with 
radius as r 2 we expect the equilibrium time at 150 au to be 



75 au, has settled and scale height enhancements are seen 



of order 4500 years. By the last frame (Fig. 13(f) ) the ma- 
jority of the disc has reached a steady state in hydrostatic 
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Figure 14. Surface density as a function of radius for the model 
disc. The solid line is the vertically isothermal disc prior to the ra- 
diative transfer code being switched on. The dotted line in shows 
the disc once is has reached a state of hydrostatic and radiative 
equilibrium. 
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Figure 15. Scale height of the model disc as a function of radius. 
The solid line is the vertically isothermal disc prior to the radia- 
tive transfer code being switched on. The dotted line in shows 
the disc once is has reached a state of hydrostatic and radiative 
equilibrium. 



and radiative equilibrium, indicating that it is not subject 
to radiative-hydrodynamic instabilities. 

Transient oscillations, as the disc adjusts towards equi- 
librium, are observed whether radiative transfer is included 
or not. They are vertical motions which decay more rapidly 
at smaller radii and hence give the appearance of propa- 
gating outwards, even though the motion is vertical. Other 
disc models using radiation hydrodynamics have found os- 
cillatory solutions with inwards propagating waves (e.g.|Min 



et al.||2009 ) which are fundamentally different to the tran- 
sient effects seen in this case. However it should be noted 



that Min et al. (2009) only found an oscillatory solution in 



their most optically thick disc, which is significantly more 
massive than our model disc, and may represent a different 
regime to that considered here. We note that the disc retains 
its azimuthal symmetry throughout the simulation and we 
find no evidence for warping instabilities. 

The surface density of the disc is shown in Fig. [14] the 
scale height is shown in Fig. ^] and the radial temperature 
profile is shown in Fig. |16| The solid lines are from the ver- 



tically isothermal disc shown in Fig. 13(a) and the dotted 



lines are the disc at the final time step shown in Fig. 13(f) 



The scale height of the disc was determined by allocating 
particles to 100 radial bins and fitting a Gaussian to the 
vertical density profile; the scale height is taken to be the 
standard deviation of the fitted Gaussian. The radial tem- 
perature profile is the average temperature within one scale 
height of the mid-plane for each radial bin. The majority 
of the surface density profile is well represented by a power 
law and does not change significantly as a result of includ- 
ing radiative transfer. Figures [l5| and [T6] show that the disc 
is hotter close to the mid-plane, except within a region be- 
tween 10-20 au of the centre, but has a smaller scale height. 
Figure ^3] shows that the vertical temperature gradient in 
the disc is such that temperature increases with height out 
of the mid-plane. Although the disc is hotter when radiative 
transfer is included the vertical temperature gradient tends 
to reduce the pressure gradient which is required for hydro- 




Figure 16. Average temperature within one scale height of the 
mid plane as a function of radius for the model disc. The solid 
line is the vertically isothermal disc prior to the radiative transfer 
code being switched on. The dotted line in shows the disc once is 
has reached a state of hydrostatic and radiative equilibrium. 



static support resulting in a smaller scale height. A slight 
drop in scale height is evident, just outwards from the disc 
inner edge, which is indicative of a puffed-up inner edge. 
Dullemond et al. (2001) find such a puffed-up inner edge 



when the disc is truncated by dust sublimation. Although 
the inner edge in our model is not due to dust sublimation, 
and is found at a greater radius than the Dullem ond et al.j 



(2001 ) case, we are still seeing a puffed up inner edge due to 
irradiation of an optically thick inner rim. 

A notable feature in the surface density plot is the drop 
in surface density at small radii. This is due to the accretion 
of material by the central sink particle and will tend to move 
the optically thick inner edge out to larger radii. Figure 
shows the SPH density in the central region of the disc at the 
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SPH particles the modelled SEDs show significant devia- 
tions from the benchmark results. 

The combined SPH-radiative transfer code was able to 
evolve a circumstellar disc to a state of radiative and hydro- 
static equilibrium. As the disc settles towards an equilibrium 
state, vertical motions are seen which give the appearance 
of an outward propagating scale height fluctuation. This is 
a vertical motion which decays more rapidly at smaller radii 
and is seen whether radiative transfer is included or not. 
These are not the same as the instabilities seen in other ra- 
diation hydrodynamics disc models (Min et al. 2009). One 
consequence of these fluctuations is that they can intercept 
radiation from the central star which can then be verti- 
cally transported into the mid-plane of the disc hence even 
thought the fluctuations are of a small amplitude they exert 
an influence on the mid-plane temperature. 



Figure 17. Kernel integrated plot of SPH density in the central 
region of the disc. This plot is from the last time step of the 
simulation and shows that the inner edge of the disc has moved 
outwards from its initial location at 1 au due to the accretion of 
material by the sink particle. 



last time step of the simulation. The inner edge of the disc 
has moved outwards from its initial value of 1 au. For the 
purposes of this model we still have an optically thick edge 
which will shadow the disc but the use of an accreting sink 
particle is likely to be inappropriate if the exact location of 
the disc inner edge needs to remain fixed. 



6 CONCLUSIONS 

We have presented results from a hybrid SPH and AMR 
grid-based radiative transfer code, which has been used to 
model a non-self gravitating circumstellar disc. Validation 
tests of the method have revealed a number of factors which 
must be taken into consideration when setting up a model 
of this type. 

Although we have found a method for constructing an 
AMR grid from SPH particles which is robust, in terms of 
the total mass on the grid, we found that it was necessary 
to add additional grid refinement in the central region of 
the disc to allow the inner edge to be well represented. The 
inner edge of the disc has a large effect on the temperature 
distribution in the disc and on the emergent SEDs. Conse- 
quently the accuracy of the solution is dependent upon the 
accurate representation of a region which constitutes a small 
fraction of the disc by volume. As the method combines 
both SPH and grid elements it is important that both these 
components have sufficient resolution to represent the disc 
inner edge. This presents a particular challenge for a SPH 
code which uses equal mass particles but is achievable with 
a large number of particles. Accurately modelling the SED 
using this method is considerably more demanding than ac- 
curately representing the temperature structure within the 
disc. Temperature errors can be reduced to a manageable 
level using a realistic number of SPH particles, yielding a 
temperature distribution suitable for use in a hydrodynam- 
ics calculation, however even with a disc composed of 10 8 
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